/*!
 **************************************************************************
 * \file convdiff2.cpp
 * This code simulate Natural Convection in a unit cube.
 *
 ***************************************************************************
 * \author    Luis M. de la Cruz [ Sat May 23 12:06:36 CDT 2009 ]
 ***************************************************************************
 * The equations to be solved are of the form:
 * \f[
 * \frac{\partial \phi}{\partial t} + \frac{\partial}{\partial x_j} 
 * \Big(u_j\phi\Big) =
 * \frac{\partial}{\partial x_j}\left(\Gamma\frac{\partial\phi}{\partial x_j}\right) + S,
 * \,\,\,\,\,\,\,\,\textrm{for  } j = 1,2,3,
 * \f]
 *
 * where \f$ \phi \f$ is a scalar variable \f$T, \rho, u_1, u_2\f$ or \f$u_3 \f$.
 * \f$ \Gamma \f$ is a diffusion coefficient and \f$ S \f$ is the source term.
 *
 * In this case we have an energy equation in terms of temperature \f$T\f$, 
 * coupled with the Navier-Stokes equations. 
 * The coupled equations are solved with the SIMPLEC strategy.
 * 
 * These equation are solved in a unit cube 
 * \f$ x, y, z \in [0,1] \times [0,1] \times [0,1] \f$ as shown in
 * the next figure:
 * \image html convection3.png "Domain in 3D" width=10cm
 * \image latex convection3.eps "Domain in 3D" width=10cm
 * Two opposite walls have constant temperatures but different. The other walls
 * are adiabatic. The no-slip condition is applied on all walls.
 *
 *  To compile and run this example type the next commands:
 * \par
   \verbatim
   % make
   % ./convdiff2 \endverbatim
 * 
 ***************************************************************************
 */

#include <iostream>
#include <cmath>

#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/Quick_CoDi.hpp"
#include "Equations/Momentum_XCoDi.hpp"
#include "Schemes/Quick_XCoDi.hpp"
#include "Equations/Momentum_YCoDi.hpp"
#include "Schemes/Quick_YCoDi.hpp"
#include "Equations/Momentum_ZCoDi.hpp"
#include "Schemes/Quick_ZCoDi.hpp"
#include "Equations/PressureCorrection.hpp"
#include "Schemes/Simplec.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;

// Scalar field in one dimension
typedef TunaArray<double,3>::huge ScalarField3D;

template<class T_e, class T_x, class T_y, class T_z, class T_p>
double SIMPLEC(T_e &, T_x &, T_y &, T_z &, T_p &, int, double);

double errorT, errorX, errorY, errorZ;

int main() 
{
    double left_wall, right_wall;
    double length_x, length_y, length_z;
    int num_nodes_x, num_nodes_y, num_nodes_z;
    double tolerance, dt;    
    int max_iter,max_time_steps;
    double Prandtl, Rayleigh;    
    int frequency;

// Reading data from the input file
    std::ifstream input_cin ("input2");
    input_cin >> length_x
	      >> length_y
	      >> length_z
	      >> num_nodes_x
	      >> num_nodes_y
	      >> num_nodes_z
	      >> left_wall
	      >> right_wall 
	      >> dt
	      >> max_time_steps
	      >> tolerance
	      >> max_iter
	      >> Prandtl
	      >> Rayleigh
	      >> frequency;


    SparseMatrix< Diagonal <double, 3> > A(num_nodes_x, num_nodes_y, num_nodes_z); 

    ScalarField3D             b(num_nodes_x, num_nodes_y, num_nodes_z); 

    StructuredMesh<Uniform< double, 3> > mesh(length_x, num_nodes_x,
					      length_y, num_nodes_y,
					      length_z, num_nodes_z);

    double dx = mesh.getDelta(X);
    double dy = mesh.getDelta(Y);
    double dz = mesh.getDelta(Z);
    int num_vols_x = mesh.getExtentVolumes(X);
    int num_vols_y = mesh.getExtentVolumes(Y);  
    int num_vols_z = mesh.getExtentVolumes(Z);  

    ScalarField3D T(mesh.getExtentVolumes());
    ScalarField3D p(mesh.getExtentVolumes());
    ScalarField3D pp(mesh.getExtentVolumes());
    ScalarField3D u(mesh.getExtentVolumes());      // u-velocity
    ScalarField3D v(mesh.getExtentVolumes());      // v-velocity
    ScalarField3D w(mesh.getExtentVolumes());      // w-velocity
    ScalarField3D us(num_nodes_x, num_vols_y , num_vols_z); // u staggered
    ScalarField3D vs(num_vols_x , num_nodes_y, num_vols_z); // v staggered
    ScalarField3D ws(num_vols_x , num_vols_y , num_nodes_z); // w staggered

    ScalarField3D Tn(mesh.getExtentNodes());
    ScalarField3D pn(mesh.getExtentNodes());
    ScalarField3D un(mesh.getExtentNodes());
    ScalarField3D vn(mesh.getExtentNodes());
    ScalarField3D wn(mesh.getExtentNodes());
//
// Valor de los nodes en las fronteras.
//
    Range all = Range::all();
    T(T.lbound(firstDim), all, all) = left_wall;   // Left
    T(T.ubound(firstDim), all, all) = right_wall;    // Rigth

    ScalarEquation<Quick_CoDi<double,3> > energy(T, A, b, mesh.getDeltas());
    energy.setDeltaTime(dt);
    energy.setNeumann(TOP_WALL);
    energy.setNeumann(BOTTOM_WALL);
    energy.setDirichlet(LEFT_WALL, left_wall);
    energy.setDirichlet(RIGHT_WALL, right_wall);
    energy.setNeumann(FRONT_WALL);
    energy.setNeumann(BACK_WALL);    
    energy.setUvelocity(us);
    energy.setVvelocity(vs);
    energy.setWvelocity(ws);
    energy.print();

    Momentum_XCoDi<Quick_XCoDi<double, 3> > mom_x(us, A, b, mesh.getDeltas());
    mom_x.setGamma(Prandtl);
    mom_x.setDeltaTime(dt);
    mom_x.setDirichlet(LEFT_WALL);
    mom_x.setDirichlet(RIGHT_WALL);
    mom_x.setDirichlet(TOP_WALL);
    mom_x.setDirichlet(BOTTOM_WALL);
    mom_x.setDirichlet(FRONT_WALL);
    mom_x.setDirichlet(BACK_WALL);
    mom_x.setUvelocity(us);
    mom_x.setVvelocity(vs);
    mom_x.setWvelocity(ws);
    mom_x.setPressure(p);
    mom_x.print();
    
    Momentum_YCoDi<Quick_YCoDi<double, 3> > mom_y(vs, A, b, mesh.getDeltas());
    mom_y.setGamma(Prandtl);
    mom_y.setDeltaTime(dt);
    mom_y.setDirichlet(LEFT_WALL);
    mom_y.setDirichlet(RIGHT_WALL);
    mom_y.setDirichlet(TOP_WALL);
    mom_y.setDirichlet(BOTTOM_WALL);
    mom_y.setDirichlet(FRONT_WALL);
    mom_y.setDirichlet(BACK_WALL);
    mom_y.setUvelocity(us);
    mom_y.setVvelocity(vs);
    mom_y.setWvelocity(ws);
    mom_y.setPressure(p);
    mom_y.setTemperature(T);
    mom_y.setRayleigh(Rayleigh);
    mom_y.print();

    Momentum_ZCoDi<Quick_ZCoDi<double, 3> > mom_z(ws, A, b, mesh.getDeltas());
    mom_z.setGamma(Prandtl);
    mom_z.setDeltaTime(dt);
    mom_z.setDirichlet(LEFT_WALL);
    mom_z.setDirichlet(RIGHT_WALL);
    mom_z.setDirichlet(TOP_WALL);
    mom_z.setDirichlet(BOTTOM_WALL);
    mom_z.setDirichlet(FRONT_WALL);
    mom_z.setDirichlet(BACK_WALL);
    mom_z.setUvelocity(us);
    mom_z.setVvelocity(vs);
    mom_z.setWvelocity(ws);
    mom_z.setPressure(p);
    mom_z.print();

    PressureCorrection<Simplec<double, 3> > press(pp, A, b, mesh.getDeltas());
    press.applyBounds(1, num_nodes_x-1, 1, num_nodes_y-1, 1, num_nodes_z-1);
    press.setUvelocity(us);
    press.setVvelocity(vs);
    press.setWvelocity(ws);
    press.set_du(mom_x.get_du());
    press.set_dv(mom_y.get_dv());
    press.set_dw(mom_z.get_dw());
    press.setPressure(p);
    press.print();    
    
    ofstream sorsumfile("./Data2/sorsum"), 
      nusseltH("./Data2/nuh"), nusseltC("./Data2/nuc"),
      errorTfile("./Data2/errorT"), errorXfile("./Data2/errorX"), 
      errorYfile("./Data2/errorY"), errorZfile("./Data2/errorZ");
    
    double sorsum, nusselt, nuh, nuc;
    int iteration = 0;
    
    NumUtils::interpolateToNodes(Tn, T);
    NumUtils::interpolateToNodes(pn, p);
    NumUtils::interpolateToNodesUs(un, us);
    NumUtils::interpolateToNodesVs(vn, vs);
    NumUtils::interpolateToNodesVs(wn, ws);
    /***
    InOut::writeToFile_DX(Tn, iteration, "./Data2/temp.", dx, dy); 
    InOut::writeToFile_DX(pn, iteration, "./Data2/pres.", dx, dy); 
    InOut::writeToFile_DX(un, vn, wn, iteration, "./Data2/velc.", dx, dy);
    ***/

    for(iteration = 1; iteration <= max_time_steps; ++iteration) {

        sorsum = SIMPLEC(energy, mom_x, mom_y, mom_z, press, 
			 max_iter, tolerance);


	if( iteration % frequency == 0 ) {
	    NumUtils::interpolateToNodes(Tn, T);
	    NumUtils::interpolateToNodes(pn, p);
	    NumUtils::interpolateToNodesUs(un, us);
	    NumUtils::interpolateToNodesVs(vn, vs);	    
	    NumUtils::interpolateToNodesVs(wn, ws);	    
	    // InOut::writeToFile_AM(Tn, iteration, "./Data2/temp.", dx, dy, dz); 
	    InOut::writeToFile_DX(Tn, iteration, "./Data2/temp.", dx, dy, dz); 
	    InOut::writeToFile_DX(pn, iteration, "./Data2/pres.", dx, dy, dz); 
	    InOut::writeToFile_DX(un, vn, wn, iteration, "./Data2/velc.", dx, dy, dz);
	}

	nuh = NumUtils::calcNusseltVert(T, dx/2.0, dy, dz, 0);
	nuc = NumUtils::calcNusseltVert(T, dx/2.0, dy, dz, num_nodes_x-1);
	sorsumfile << iteration << "\t" << sorsum << endl;
	errorTfile << iteration << "\t" << errorT << endl;
	errorXfile << iteration << "\t" << errorX << endl;
	errorYfile << iteration << "\t" << errorY << endl;
	errorZfile << iteration << "\t" << errorZ << endl;
	nusseltH << iteration << "\t" << nuh << endl;
	nusseltC << iteration << "\t" << nuc << endl;
	
	cout << "=> Iter = " << iteration 
	     << " s = " << sorsum
	     << " e_T = " << errorT
	     << " e_X = " << errorX
	     << " e_Y = " << errorY
	     << " e_Z = " << errorZ
	     << " nuh = " << nuh
	     << " nuc = " << nuc << endl;
    }
   
    std::cout << "\n\n NOTE: The data were written in the [./Data2] directory in DX format";

    return 0;
}


template<class T_e, class T_x, class T_y, class T_z,
	 class T_p>
double SIMPLEC(T_e &energy, T_x &mom_x, T_y &mom_y, T_z &mom_z, T_p &press, 
	       int max_iter, double tolerance) 
{
    double sorsum = 10.0, tol_simplec = 1e-02;
    int counter = 0;

    int pausa;

    while ( (sorsum > tol_simplec) && (counter < 20) ) {

	energy.calcCoefficients();
	Solver::TDMA3D(energy, tolerance, max_iter);
	errorT = energy.calcErrorL2();
	energy.update();

	mom_x.calcCoefficients();
	Solver::TDMA3D(mom_x, tolerance, max_iter);
	errorX = mom_x.calcErrorL2();
	mom_x.update();

	mom_y.calcCoefficients();
	Solver::TDMA3D(mom_y, tolerance, max_iter);
	errorY = mom_y.calcErrorL2();
	mom_y.update();

	mom_z.calcCoefficients();
	Solver::TDMA3D(mom_z, tolerance, max_iter);
	errorZ = mom_z.calcErrorL2();
	mom_z.update();

	press.calcCoefficients();
	Solver::TDMA3D(press, tolerance, max_iter);
	press.correction();	
	sorsum = fabs( press.calcSorsum() );

	++counter;
    }

    return sorsum;
    
}





